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1.  INTRODUCTION 

The  computer  code  A32CODE  performs  automatic  analysis  of  the  AFRPL 
single-color  transmlssooeter  data  to  obtain  the  volume-to-surface  mean  radius 
a32  of  plu®E  particles  and  estimates  of  Its  error  caused  by  uncertainty  In 
size  distribution,  complex  index  of  refraction,  and  experimental  data*  The 
code  consists  of  two  parts:  the  actual  program  (A32CODE)  and  subprograms  for 
analyzing  the  data,  and  files  of  Mie  scattering  parameters  (UQFILE  and 
BQFILE). 

The  method  of  analysis  has  been  presented  in  detail  in  Section  4.1  of 
Ref.  1  and  is  only  briefly  reviewed  here  in  the  discussions  of  how  the  code 
works.  The  input  data  required  by  the  code  are  defined  in  Table  1.  With  the 
data  of  the  first  card  and  the  data  of  either  UQFILE  (if  DISTR  “  U)  or  BQFILE 
(if  DISTR  ■  B),  the  upper  and  lower  error  bounds  of  Q/a^  versus  are 
constructed  and  stored.  With  the  data  from  the  second  and  subsequent  cards, 
the  experimentally  determined  value  of  Q/a^  is  computed  and  used  with  the 
constructed  curves  to  determine  aj2  and  its  error.  An  example  application  is 
described  in  a  later  section. 
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Table  1.  A32CODE  Card  Input  Format 


loading  and  error  (g/cm  ) 


2.  MIE  PARAMETER  PILES 


The  two  flies  UQPILE  and  BQPILE  contain  the  Mle  scattering 
parameter  Q/x^  averaged  over  all  possible  unlaodal  and  bleodal  rectangular 
size  distributions,  respectively,  for  which  x32  is  fixed  at  the  tabulated 
value.  The  files  contain  entries  for  the  54  pairs  of  n,  k  values  laplled  by 

n  -  1.60,  1.65,  1.70,  1.75,  1.80,  1.85,  1.90,  1.95,  2.00 
•c  -  0,  2  x  10"4,  2  x  10-3,  2  x  10"2,  0.2,  0.5 

The  structure  of  ei^^r  file  is  defined  by  the  FORTRAN  read  sequence  of 
subroutine  READQ  (a  listing  of  the  code  is  given  in  the  Appendix).  NINDEX(N) 
(N-1,9)  and  KINDEX(K)  (K-l ,6)  are  the  arrays  just  listed.  X32(I)  (1-1,41)  is 
the  array 


x32(i)  -  100,1x32(i-l);  1-1,41}  x32(l)  -  0.1 

QBAR(N,K,1)  is  the  value  of  Q/x32  averaged  over  the  size  distribution, 
QSIG(N,K,1)  is  the  standard  deviation  of  variation  about  the  mean  value, 
QMIN(N,K,I)  is  the  oiniaua  value  in  the  variation,  and  (J1AX(N,K,I)  is  the 
eaxiBua  value  in  the  variation. 


3.  A32CODB  DESCRIPTION 


The  first  step  In  A32CODE  Is  to  read  the  input  date  of  the  first  data 
card  and  the  data  of  either  UQPILE  or  BQPILE.  The  latter  is  done  by 
subroutine  READQ.  The  data  of  the  first  card  specify  the  degree  of  ignorance 
assured  in  n  and  x,  the  wavelength,  and  whether  the  analysis  is  to  be  done 
using  a  unlaodal  or  biaodal  rectangular  particle  size  distribution  (see  Table 
1  for  details)* 

The  next  step  is  to  effect  a  variation  over  n  and  x.  The  ranges  of  n 
and  k  are  divided  into  NN-1  and  KK-1  equal  size  subintervals  and  the  five  Q~ 
var tables  described  above  are  computed  at  the  NN*  KK  grid  points  for  all  41 
values  of  X32  (ttie  code  uses  NN»KK“5,  but  these  nunbers  can  easily  be 
changed)*  The  computation  of  the  Q  values  consists  of  interpolations  on  the 
values  in  either  the  UQPILE  or  BQPILE.  The  interpolation  is  of  the  fora 

Log  Q  -  (A-fBn)  +  (C+On)x  (1) 

and  is  performed  in  subroutine  QINT*  At  each  of  the  41  values  of  X32,  two 
upper  and  two  lower  bounds  of  Q/«32  are  determined.  These  are 

CHAX1  -  max  (QBAR+QSIG)  c 

QMAX2  -  max  (QMAX)  c  (2) 

QMIN1  -  min  (QBAR-Q6IG)  c 
QMIN2  -  min  (QMIN)  c 

where  max  (  )  and  min  (  )  denote  the  maximum  and  minimum  values,  respectively. 
In  the  NNx  KK  set  of  variation  values.  QNAX1  and  QHIN1  represent  ’’optimistic" 
error  bounds  in  that  they  are  the  traditional  root-mean-square  (rv) 
deviations  from  the  mean.  QMAX2  and  QMIN2  represent  "pessimistic"  error 
bounds  in  that  they  are  the  absolute  maximum  and  minimum  values  in  the 
variation  over  index  of  refraction  and  size  distribution.  If  QMIN1  is  less 
than  or  equal  to  zero,  it  is  replaced  with  CPIN2.  The  constant  c  is  2w/X  and 
transforms  Q/x^  *nto  ^*32  un*c 

The  end  result  of  the  variation  then,  is  two  sets  of  upper  and  lower 
bound  curves  versus  832*  In  the  code,  these  curves  are  denoted  by 
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BL(m,i)  I 

i-1,41 

BU(m,i) J 

where  L  implies  "lower",  (J  lap  lies  "upper",  ur*l  Implies  ras  error  sad  op2 
Implies  full-range  error.  These  curves  are  essentially  the  bound  curves  of, 
for  example.  Fig.  31  of  Ref.  1.  The  four  curves  are  plotted  versus  a^  In 
subroutine  PLQTB  (note,  X32  “  )• 

For  subsequent  analysis,  it  is  convenient  to  have  an  algebraic 
representation  of  these  curves.  Here,  a  quadratic  fit  is  made  on  the 
logarithm  of  the  curves.  In  each  Interval  a^  to  a£+i  with  1-1,40,  the  curves 
are  approximated  by  (B(x)  is  any  one  of  the  bound  curves  and  a  is  832) 


where 


for  i»l,39,  and 


Bi(a)  -  e 


+  B^n  (a/aj)  +  a^ln  (a/ai) 


(3) 


1 


1  21n2(a1+1/ai) 


in 

»i+l 
_  4 


*i  21n 


l+r  i 

T,  -  In  B 


V  Bi+2 


(4a) 

(4b) 

(4c) 


“l  '  °39  81  '  639  \  m  hi 


(5) 


for  i-40. 


In  the  code,  these  coefficient  arrays  are  denoted  by 


ACL(a,i) 
BCL(m,i) 
CCL(a.i) 
ACU(n,i) 
BCU(n,l) 
CCU(m.i)  , 


i-1,40 

m-1,2  (i.e.,  rms  or  full-range  error) 
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This  curve  fitting  is  done  in  subroutine  CURFIT.  In  addition  to  the  curve 
fitting,  the  peak  value  of  each  curve  is  computed  in  CURFIT  and  stored  as 

YLMAX(m) 

YUMAX(a) 

At  this  point  in  the  code,  a  loop  is  initiated  over  experimental  input. 
One  card  at  a  time  is  read  in,  the  data  contained  on  it  are  analyzed,  and  the 
results  for  mean  particle  size  and  error  printed  out.  The  loop  continues 
until  an  end-of-information  card  is  reached,  at  which  point,  the  code  stops. 

The  experimental  input  data  required  are  described  in  Table  1.  From 
these  data,  the  experimental  value  of  Q/a^  (denoted  by  Y  in  the  code)  is 
computed  from  (see  Eq.  (42)  in  Ref.  1] 

Y(in  pm"1)  -  -  -3^7-  lO"*  In  r  (6) 

m 

Its  experimental  error  is  determined  from 


The  experimental  bounds  on  Q/a^  are  then  defined  by 

YP  -  Y  +  DY  (8) 

YM  -  Y  -  DY 

At  this  point,  an  analysis  is  made  of  where  these  two  bounds  (le.,  YP  and 
YM)  intersect  the  BO  and  BL  curves.  (A  loop  is  made  over  m-1  and  2  so  that 
the  following  is  done  for  both  the  rms  and  full-range  error  bound  cases).  The 
primary  differentiation  on  the  manner  of  analysis  is  made  by  testing  the 
magnitudes  of  YP  and  YM  with  YUMAX  and  YIMAX.  The  most  general  case  occurs 
when  YP  is  less  than  YLMAX,  that  is,  when  the  upper  experimental  bound  lies 

below  the  peak  of  the  lower  bound  curve.  This  condition  is  shown  in  Fig.  1. 

In  this  case  there  are  two  well-defined  regions  of  the  a^  axis  where  the 
experimental  data  are  consistent  with  the  analysis  curves.  They  lie  on  either 
side  of  the  peaks  of  the  curves  BU  and  BL  and  are  designated  here  as  lying  in 
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Che  "upper"  or  ’’lower"  branch  of  the  analysis  curves.  The  eight  points  of 
intersection  of  YP  and  YM  with  BD  end  BL  shown  in  Fig.  1  are  computed  in 
subroutine  ROOTS  (described  later).  The  retrieved  a^  and  errors  are  defined 
in  terns  of  the  032  coordinates  of  these  points  by 


*32 


Aa  - 

32 


T  £aUl  +  aLl  +  aUl  +  *L1  1 

1  r  +  +  ~  —  1 

T  laL2  +  aU2  +  aL2+  aU2  J 
2  ^*L1  “  Vl) 

1  ^*U2  "  *L2  ^ 


on  lower  branch 
on  upper  branch 
on  lower  branch 
on  upper  branch 


(9a) 


(9b) 


(The  code  syabols  for  these  intersection  points  are  AD1P,  AL1P,  AU1M,  ALIM, 
AL2P,  AU2P,  AL2M  and  AU2M). 

In  no8t  cases,  it  will  be  known  fron  experinental  conditions  which  of  the 
two  retrieved  values  is  relevant.  For  those  cases  where  it  is  not  known,  the 
code  also  computes  a  "branch  undetermined"  retrieval  value  for  a32.  In 
principle,  this  value  and  its  error  could  be  determined  from  the  four 
intersections  of  TP  and  YM  and  BU  by 


*32 


1  r*Dl  *  *01 

2  I  2 


,  *U2  *  *P2  1 

2  J 


(10) 


**32  "  2  (*U2  "  *U1 ) 


In  practice,  however,  these  values  are  more  conveniently  computed  in  the  code 
from  the  Intersections  of  the  mean  experimental  value  T  with  the  upper  bound 
curve  BU  (Fig.  1)  by 


a32  "  2  <al  +  *2^ 
A*32  "  2  ^*2  ~  *1^ 


(ID 


(The  code  symbols  for  these  Intersection  points  are  A1  and  A2.) 


If  Che  upper  experimental  bound  YP  la  greeter  than  the  peak  YLMAX  of  the 
lower  bound  curve,  there  are  no  well "defined  "upper"  and  "lower"  branch 
regions.  For  this  case,  a  "branch  undefined"  analysis  is  made  exactly  as  for 
a  "branch  undetermined"  case.  That  is,  the  intersections  of  T  with  BO  are 
determined  and  Eqs.  (11)  used  to  define  the  retrieval  a^  and  error.  A  minor 
variation  is  made  if  T  is  greater  than  TtiMAX  (so  that  there  are  no 
intersections  of  T  with  BU)  —  the  value  of  YM  is  used  in  place  of  T.  If  YM 
is  greater  the  YOMAX,  no  analysis  is  possible.  (Also,  analysis  is  skipped  if 
YM  is  less  than  the  lowest  value  of  BL  in  the  entire  range.) 

For  all  cases  of  analysis,  the  intersection  of  a  line  Y  (l.e.,  Y,  YP,  or 
YM)  with  a  bound  curve  B  (i.e.,  BU  or  BL)  is  asde  in  subroutine  ROOTS  using 
the  algebraic  representation  for  B  defined  by  Eqs ■  (3)-(5).  In  each  of  the  40 
subintervals ,  a  test  is  made  to  see  if  there  are  any  intersections  of  Y  with 
the  quadratic  form  of  Eq.  (3).  There  will  be  no  real  intersections  if 

F  "  (Yj-ln  Y)  <  0 

and  two  Intersections  (they  may  be  equal)  otherwise.  These  intersections  are 

.<*•»  ..t .  <**  V/2ai 

A  further  check  is  made  to  see  If  these  roots  occur  in  the  Interval  aj  to 
al+l*  they  do,  they  are  saved.  After  all  roots  in  all  40  subintervals 
have  been  found  (there  may  be  many  roots  If  B  fluctuates  up  and  down  a  lot)  a 
search  is  made  for  the  smallest  and  largest.  These  two  are  returned  to  the 
main  code  as  the  relevant  intersection  points. 


4.  EXAMPLE  APPLICATION 


Example  run*  of  the  code  Mere  nude  for  the  eeee  where  Ignorance  of  the 
particulate  Index  of  refraction  le  specified  by 

1.7  <  n  <  1.8 

0  <  k  <  0.02 

and  X  le  the  AFRPL  experimental  value  X  ■  0.S14S  m. 

Both  the  unimodal  and  blaodal  alee  distribution  were  used.  Two  cases  of 
experimental  data  were  run.  In  both, 

t  •  0.880 
L  -  10  cm 

Cm  -  5.5  x  10“7g/cm3 
d  -  3.7  g/cm3 


In  the  first  case,  these  values  were  assumed  to  be  exact  and  the  Input  errors 
were  set  to  zero.  The  resulting  error  In  032  !■  then  due  soley  to  Ignorance 
of  size  distribution  and  Index  of  refraction.  In  the  second  case,  the 
following  experimental  errors  were  assigned: 

At  -  0.03  (~2Z  error  In  I  and  I  ;  t  ■  I/I  ) 

o  o 

AL  -  1  cm  (10Z) 

AC—  ■  5.5  x  10  *  g/cm^  (10Z) 

Ad  "  0.4  g/cm3  (~10Z) 

Code  results  for  these  two  cases  with  the  unlmodal  size  distribution  are 
shown  In  Tables  2  and  3  and  Figs.  2  and  3.  Table  2  lists  the  Q/a^  error 
bounds.  AVEASIG  are  the  "optimistic"  error  bounds  wherein  the  variation  over 
size  distribution  Is  measured  by  the  standard  deviation  of  Q/a^*  MAX  and  MIN 
are  the  "pessimistic"  error  bounds  wherein  the  variation  over  size 
distribution  Is  measured  by  the  total  range  of  variation  of  Q/a^*  The  choice 
of  whether  to  be  optimistic  or  pessimistic  is  a  user  option.  These  bound 
results  are  plotted  In  Figs.  2  and  3.  (Note:  the  Aerospace  plotting  routines 


Tabic  2.  Q/*32  Bounds  Listing  for  Dniaodal  Distribution 

A32COOE  RESULTS 

NLIHITS  NX-  1.700E*00  N2»  1.000E*00 

K LIMITS  Kl=  0.  K2=  2.S00E-02 

HAVE LENGTH*  .5X45  MICRON 
DISTRIBUTIONS  U 


Q/A32  BOUNDS  <UNIT=1/HU> 


X 

A32IHU) 

X32 

AVE«SI8 

AVE-SIO 

MAX 

HXN 

X 

S.1B9E-03 

1 . OOOE-Ol 

4.271E-01 

S.241E-03 

4.273E-01 

4.S7SE-03 

2 

1.032E-02 

1.260E-01 

4.3I9E-01 

1.0S0E-02 

4.361E-01 

9.770E-03 

3 

1.294E-02 

1.580E-01 

4.S21E-01 

2.074E-02 

4.526E-01 

1.930E-02 

4 

1.630E-02 

2. 000E-01 

4.646E-01 

4.222E-02 

4.854E-01 

3.925E-02 

5 

2.055E-02 

2.S10E-01 

5.454E-01 

S.3S9E-02 

5.469E-01 

7.794E-02 

6 

2.5S6E-02 

3.160E-01 

6.O47E-01 

1.6S7E-01 

6.883E-01 

1.567E-01 

7 

3.259E-02 

3.980E-01 

9.777E-01 

3.411E-01 

9.049E-01 

3.163E-01 

6 

4. 102E-02 

5.010E-01 

1.567E+00 

6.908E-01 

1.581E*00 

6.409E-01 

9 

S.167E-02 

6.310E-01 

2.747E*00 

1.402E*00 

2 . 776E*00 

1.302E*Q0 

10 

6.S02E-02 

7.940E-01 

4.955E*00 

2.795E+00 

5.005E*00 

2.616E«00 

11 

8.189E-02 

1.000E*00 

6.621E*00 

S.242E*00 

0.712E400 

5.029E*00 

12 

1.032E-01 

1.260E*00 

1.546E*01 

9.099E*00 

1.553E+01 

8.649E+00 

13 

1.294E-01 

1 . 580E*00 

2.098E*01 

1.5471*01 

2 . 231E*01 

1.545E*01 

14 

1.63SE-01 

2.000E+00 

2.258E+01 

1.025E*01 

2 . 304E+01 

1.809E*01 

IS 

2.055E-01 

2.510E*00 

2.157E*01 

1.760E*01 

2 . 262E+01 

1.768E*01 

16 

2.583E-01 

3. 160E*00 

1.709EO1 

1.311E*01 

1.889E+01 

1.296E*01 

17 

3.259E-01 

3.980E+00 

1.067E+01 

8.138E*00 

1.140E*01 

7.843E+00 

IB 

4.102E-01 

5.010E+00 

6 . 724EA00 

5.015E*00 

6.645E*00 

4. 148E+00 

19 

5.167E-01 

6 . 310E+00 

5.494(400 

4.641E+00 

5.679E+00 

4.312E+00 

20 

6.502E-Q1 

7.940E»00 

4.482E+00 

3.401E*00 

5.154E+00 

3.661E«00 

21 

0.109E-01 

1.000E+01 

3.167E+00 

2.636E*00 

3.081E*00 

2.225E*00 

22 

1.032E+00 

1 . 260E*01 

2.447E+00 

2.208E+00 

2.628E*00 

2.121E*00 

23 

1.294E+00 

1.580E*01 

1.862E*00 

1.766E*00 

1.930E+00 

1 . 744E*00 

24 

1.638E*00 

2.000E+01 

1.415E+00 

1.354E+00 

1.450E+00 

1 . 305E*00 

25 

2 .055E*00 

2.510E+01 

1.132E+00 

1.071E*00 

1.206E+00 

1.050E*00 

26 

2.588E*00 

3.160E+01 

S.564E-01 

8.357E-01 

6.602E-01 

S.147E-01 

27 

3.259E+00 

3.980E+01 

6.746E-01 

6.573E-01 

6.079E-01 

6.443E-01 

2« 

4.102E*00 

5.010E+01 

5.312E-01 

5.154E-01 

5.433E-01 

5.019E-01 

29 

5.167E*00 

6 . 310E+01 

4.136E-01 

4.076E-01 

4.164E-01 

4.015E-01 

30 

6.502E*C0 

7.940E*01 

3. 286E-01 

3.230E-01 

3.351E-01 

3. 209E-01 

31 

8. 189E+00 

1.000E*02 

2.576E-01 

2.53OE-01 

2.606E-01 

2.510E-01 

32 

1.032E401 

1 . 260E+02 

2.026E-01 

2.I04E-01 

2.031E-01 

1.9S4E-01 

33 

1 . 294E  *01 

1 . 580E*02 

1.608E-01 

1.591E-01 

1.621E-01 

1.577E-01 

34 

1.638E+01 

2 . 000E*02 

1.265E-01 

1.255E-01 

1.273E-01 

1 . 255E-01 

35 

2.055E*01 

2.510E*02 

9.999E-02 

9.958E-02 

1.002E-01 

9. 930E-02 

36 

2.58SE*01 

3. 160E*02 

7.907E-02 

7.883E-02 

7.915E-02 

7.865E-02 

37 

3. 259E*01 

3. 980E*02 

6 . 269E-02 

6.242E-02 

6 . 281E-02 

6 . 226E-02 

30 

4. 102E*01 

5.C10E+02 

4.959E-02 

4.951E-02 

4.968E-02 

4.951E-02 

39 

5. 167E*01 

6 . 310E+02 

3.929E-02 

3.922E-02 

3.932E-02 

3. 921E-02 

40 

6.502E+01 

7.940E*02 

3.116E-02 

3.110E-02 

3.114E-02 

3.106E-02 

41 

8. 189E*01 

1.000E+03 

2.470E-02 

2.467E-02 

2.473E-02 

2.467E-02 

18 


Table  3.  Retrieval  Results  for  Unlaodal  Distribution 


A32  AMD  TOTAL  ERROR  RESULTS 

FIRST  PCW=VALUE  LB -LOWER  BRANCH 

SECOND  RCW-RMS  ERROR  UB=UPPER  BRANCH 

THIRD  ROW-PERCENT  ERROR  0B=UNDETERt1INED  BRANCH 

ZERO=NOT  APPLICABLE 

TRANS  LCCM)  CC6/CH3)  0CS/CH3)  R/A(MU-1»  A32(HU>  RMS  ERROR  BOUNDS  A32(HU>  FULL  ERROR  BOUOS 

LB  UB  OB  LB  UB  OB 


6.600E-01  1 .003E*01  S.500E-07  3.700E*00  1.147E+01  1.01E-01  2.9SE-0X  2.02E-01  1.02E-0I  3.00E-01  2.08E-01 
0  o  0.  0.  0.  1.08E-0Z  1.9ZE-02  1.12E-01  1.14E-02  2.48E-02  1.17E-01 
0;  0.  0.  0.  0.  1.07E*01  6.49E+00  5.53E+01  1.12E+01  8.25E+00  5.65E+01 


8.800E-01  1.000E+01  5.S00E-07  3.700E«00  1.147E+01 
3.000E-02  1.000E*00  5.500E-08  4.000E-01  3.676E+00 
2.667E+01  1.000E+01  1.000E+01  1.081E+01  3.2C6E+01 


1.01E-01  3.04E-01  2.02E-01  1.02E-01  3.06E-01  2.08E-01 
2.48E-02  7.00E-02  1.12E-01  2.51E-02  7.08E-02  1.17E-01 
2.45E+01  2.30E+01  5.53E*01  2.46E+01  2.31E*01  5.65E+01 
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are  lnhouae  and  not  transferable  to  other  computer  facilities;  the  plot 
routine  in  the  delivered  version  of  the  code  gets  everything  ready  to  plot, 
but  does  not  actually  plot  then, ) 

Table  3  lists  the  actual  retrieved  a^  *nd  its  error  for  both  the 
optimistic  and  pesslalstlc  error  bounds*  For  each  of  these,  three  values  are 
determined  depending  on  which  branch  (if  either)  of  the  Q/a^  curve  is 
chosen.  For  the  AFRPL  conditions,  it  is  probably  safe  to  assume  the  upper 
branch* 

Tables  4  and  S  and  Figs.  4  and  5  repeat  the  results  for  the  biaodal  else 


distribution* 


• 

Table  4. 

Q/a^2  Bounds  Listing 

for  Blaodal  Distribution 

• 

A32COOE  RESULTS 

■  "  !;■  ’ 

NLXMXTS 

Nl*  1 

.  700E+00  N2*  1.800E+00 

KLIMTS 

Kl*  0 

K2= 

2.000E-02 

NAVE LENGTH* 

.5145  MICRON 

• 

DISTRIBUTION* 

B 

Q/A32  BOUNDS  (UNIT* 1/MU) 

■ ;; 

I 

A32(MU) 

X32 

AVE*SI6 

AVE-SX0 

MAX 

MIN 

*  -'-—9 

1 

8.189E-03 

1.000E-01 

4.484E-01 

4.87SE-03 

5.659E-01 

4.875E-03 

• 

2 

1.032E-02 

1.260E-01 

4.754E-01 

9.771E-03 

7.3476-01 

9.771E-03 

3 

1 . 294E-02 

1.580E-01 

5.2856-01 

1.9106-02 

1.074E*00 

1.9106-02 

4 

1.638E-02 

2.000E-01 

6.5356-01 

3.9256-02 

1.7436*00 

3.925E-02 

5 

2.055E-02 

2.510E-01 

8.9366-01 

7.732E-02 

2.8706*00 

7.732E-02 

» 

6 

2.588E-02 

3.160E-01 

1.331E*00 

1.567E-01 

4.901E*00 

1.567E-01 

7 

3.259E-02 

3.980E-01 

2.0206*00 

3.164E-01 

7.9736*00 

3.1646-01 

8 

4.102E-02 

5.010E-01 

2.960E*00 

6.408E-01 

1.3756*01 

6.408E-01 

9 

5.1676-02 

6.310E-01 

3.9396*00 

1.3026*00 

1.291E*01 

1.302E*00 

w 

10 

6.502E-02 

7.940E-01 

5.7S9E*00 

2.6136*00 

1.4076*01 

2.613E»00 

11 

8.189E-02 

1.0006*00 

9.3056*00 

5.0746*00 

1.336E*01 

5.029E*00 

12 

1.032E-01 

1.2606*00 

1.574E*01 

8.6856*00 

1.8506*01 

5.0036*00 

13 

1.294E-01 

1.5£0E*00 

2.229E*01 

1.381E«01 

2.231E*01 

3.9096*00 

14 

1.638E-01 

2.0006*00 

2.277E+01 

1.6556*01 

2.3156*01 

5.089E*00 

15 

2.C55E-01 

2.5106*00 

2.217E*01 

1.6716*01 

2.282E*01 

6.8616*00 

16 

2.588E-01 

3.1606*00 

1.788E*01 

1.6736*01 

1.8886*01 

7.078E*00 

17 

3.259E-01 

3.980E«00 

1.114E+01 

7.6166*00 

1.430E*01 

6. 334E*00 

t  .  . 

18 

4.102E-01 

5.010E+00 

6.7216*00 

4.6096*00 

9.4696*00 

4.1526*00 

19 

5.167E-01 

6.310E*00 

5.6906*00 

4.4916*00 

7.375E*00 

3.4306*00 

20 

6.592E-01 

7.9406*00 

4.7396*00 

3.7256*00 

5.1516*00 

3.097E*00 

21 

8.189E-01 

1.000E*01 

3.1206*00 

2.4786*00 

3.905E*00 

2.2246*00 

22 

1.0326*00 

1.2606*01 

2.544E«00 

2.1706*00 

2.6946*00 

1.829E*00 

23 

1.294E*00 

1.5806*01 

1.9036*00 

1.7546*00 

2.055E*00 

1.6566*00 

24 

1.638E*00 

2.000E*01 

1.440E*00 

1.3316*00 

1.5036*00 

1.3036*00 

25 

2.055E*00 

2.510E+01 

1.163E*00 

1.0576*00 

1.206E*00 

1.039E*00 

t 

26 

2.5886*00 

3.160E*01 

8.642E-01 

8.2456-01 

9.416E-01 

8.1006-01 

27 

3.2596*00 

3.9806*01 

6.8046-01 

6.5276-01 

6.994E-01 

6.440E-01 

28 

4.102E*00 

5.0106*01 

5.356E-01 

5.0996-01 

5.4326-01 

5.016E-01 

29 

S.167E*00 

6.3106*01 

4. 151E-01 

4.0506-01 

4.2776-01 

3.9906-01 

30 

6.502E+00 

7.940E+01 

3.313E-01 

3.2206-01 

3.352E-01 

3.190E-01 

31 

8.1896*00 

1.0306*02 

2.5916-01 

2.526E-01 

2.642E-01 

2.5076-01 

32 

1.0326*01 

1.2606*02 

2.027E-01 

1.9956-01 

2.0636-01 

1.9836-01 

33 

1.294E*01 

1.560E*02 

1.6166-01 

1.5846-01 

1.6246-01 

1.5746-01 

• 

34 

1.633E*01 

2.000E+02 

1.267E-01 

1.2566-01 

1.2766-01 

1.246E-01 

35 

2.0556*01 

2.5106*02 

1.0026-01 

9.9636-02 

1.011E-01 

9.924E-02 

36 

2.588E*01 

3.160E*02 

7.915E-02 

7.8796-02 

7.9956-02 

7.856E-02 

37 

3.259E*01 

3.9806*02 

6.274E-02 

6.2366-02 

6.3346-02 

4.2206-02 

'  '  * 

38 

4.102E*01 

5.010E*02 

4.967E-02 

4.9506-02 

4.9916-02 

4.9416-02 

*  *  * 

39 

5.1676*01 

6.3106*02 

3.933E-02 

3.9226-02 

3.945E-02 

3.9156-02 

40 

6.502E*01 

7.940E*02 

3.117E-02 

3.1076-02 

3.1286-02 

3.0936-02 

41 

8.1896*01 

1.0006*03 

2.4606-02 

2.4416-02 

2.4776-02 

2.4406-02 

i- 

23 
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Bounds  -  Biaodal  Size  Distribution,  Full-Range  Error 


APPENDIX 


A32C0DE  Listing' 


100  PROGRAM  A32CQDE( INPUT* OUTPUT *TAPE2*TAPE3* TAPES* INPUT *TAPE6=0UTPUT) 

110  C 

120  DIMENSION  YLMAX<2> » YUMAXI2) 

130  DIMENSION  AUB< 2) »DUB<2> *EUB< 2) 

140  DIMENSION  ALB<2)*DLB<2)*ELB<2) 

150  DIMENSION  A06<2> *D0B<2> *E0B<2> 

160  REAL  N1 *N2*K1 *K2* NINDEX (10) *  KINDEX< 10 )  *  L 

1 70  C0MM0N/QDATA/X32<  41 ) »XXXX1 < 8873 > 

190  COMMON/QFUNC/QBARl 41 ) *QSIG<  41 ) *QMIN(41 ) *9MAX<  41 ) 

190  C0MM0N/B0UND/A32(41 ) *BU<2»41)*BL(2*41) 

200  C 
210  C 

220  C  READ  AND  LIST  INPUT  DATA 

230  READ(5* 100)  N1»N2.K1*K2*NL>DISTR 

240  WRITE<6*200> 

250  WRIT£<6*  201 )  N1*N2 

260  WR!TE<6»202>  K1*K2 

270  WRITE(6*203>  UL 

280  WRITE<6*204>  DISTR 

290  C  READ  QFILE  DATA 
200  CALL  READO(DISTR) 

310  C 

320  C  SET  UP  VARIATION  ARRAYS 
330  NN=5 

340  XK=5 

350  DELN=(N2-N1)/<NN-1.> 

360  DO  1  N- 1 »NN 

370  1  NINDEX<N)=N1+<N-1.)*DELN 

380  DELK=(K2-K1)/<KK-1.) 

390  DO  2  K=1.KK 

400  2  KINDEX(K)=K1+(K-1.)*DELK 

410  C  SET  UP  A32  ARRAY  AND  INITIALIZE  FOR  UPPER  AND  LONER  BOUNDS  SEARCH 

420  CC =6. 283 185308/ WL 

430  DO  3  1=1*41 

440  A32( I ) =X32< I )/CC 

450  DO  3  M=1  *2 

460  BU<M*I)=0. 

470  3  BL<M»I)=1.E99 

480  C  PERFORM  VARIATION  OVER  INDEX  OF  REFRACTION  AND  SET  UP  LOWER  AND 

490  C  UPPER  BOUNDS  SEARCH 

500  DO  4  N=1 *NN 

510  DO  4  K=1 *KK 

520  CALL  GINT(NINDEX(N) >KINDEX(K) > 

530  DO  4  1=1*41 

540  QMAX1=CC*(9BAR<I)+0SIGU>> 

550  QM1N1=CC*'.0BAR<I)-0SIGU>) 

560  QMAX2=CC*0MAX( I ) 

570  QM1N2=CC*QM!N<I> 

580  IF(QMIN1.LT,QMIN2)  QMIN1=0MIN2 
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590  IF (QMAX1 . ST . BU < 1 » I ) )  BUI 1 » I )=QHAX1 

600  IF<QM!N1.LT.BL!1»I>>  BL< 1 > I )=QMIN1 

610  IF(QMAX2.3T .BU(2t I ) )  B'J<2»  I  )=QMAX2 

620  IF(QMIN2.LT.BL<2»I>>  BL<2« I >«QMIN2 

630  4  CONTINUE 

640  C  PRINT  BOUNDS  RESULTS 

650  WRITE(6»300> 

660  WRITE(6»301 > 

670  DO  5  1  =  1 » 41 

680  5  WRITE<6*302>  I »A32< I ) »X32(I ) » BU< 1>I)»BL(1» I ) »BU<2» I > »BL<2» I > 

690  C  FIT  UPPER  AND  LOWER  BOUND  RESULTS  WITH  OUADRATIC  SEGMENTS 

700  DELX=0. 230258510 

710  CALL  CURFIT ( 0» 1 » 41 »DELX?  YLMAXI 1 ) ) 

720  CALL  CURFIT ( 1 » 1 »41 >DELX* YUNAX( 1 ) ) 

730  CALL  CURFIT(0»2»41 tDELXi YLMAXI2) ) 

740  CALL  CURFIT (l*2t41f DELX>YUMAX(2) ) 

750  C  PLOT  BOUND  RESULTS 
760  CALL  PLOTB(YUMAX) 

770  C 

780  C  DATA  PROCESSING  LOOP 
790  WRITE(6*40C> 

900  WRITE' 6»40t> 

810  WRITEI6*  402' 

820  WRITE<6»403> 

930  WRITE ( 6»  404  ) 

840  WRITE! 6. 405) 

850  C  READ  EXPERIMENTAL  DATA  AND  COMPUTE  EXPERIMENTAL  ERROR 
860  6  READ(5»101)  TR.DTR.Lr PL»C»DC»D»DD 

870  IFIEOF (5) )  7.8 

890  7  STOP 

890  8  ET=-100.*DTR/<TEtALQ6<TR>> 

900  EL- IOC. *DL/L 

910  EC-100. *DC/C 

920  ED=100.*DD.'t> 

930  Y=-4.E-4*D*ALOGITR)/<3.*L*C) 

940  EY=SQRT(ET*ET+EL*EL+EC*EC+ED*ED> 

950  DY-EYBY/100. 

960  YP-=  Y+DY 

970  '( M=Y-DY 

990  C  LDOF  OVER  TWO  CASES  OF  SIZE  DISTRIBUTION  ERROR  --  M*1  IMPLIES  RMS* 
990  C  M~2  IMPLIES  FULL-RANGE 
1000  DO  11  M  *  1 » 2 

1010  C  INITIALIZE  BRANCH  RESULTS 
1020  AUB(M)=0. 

1030  ALB<M)=0, 

1040  A0B(M>=0. 

1050  DUB(M)=0. 

1060  DLB(M)=0, 

1070  DOB<M>=0. 

1080  EUB(M)=0. 
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1090  ELB!  MHO. 

1100  E0B!M)=0. 

1110  C  UNDETERMINED(OR  UNKNOWN)  BRANCH  ANALYSIS 
1120  Z=Y 

1130  !F(Z.LE.YUNAX(M))  GO  TO  9 

1  HO  Z=YM 

1150  IF1Z.GT .YUHAX(N) )  GO  TO  10 

1160  9  CALL  ROOTS! l»Mi 41 >Z»A1»A2) 

1170  AOB!M)=< A24A1 )/2. 

1180  D0B1MH!  A2-A1  )/2. 

1190  EOB! M)S100.*D0B<N)/A0B!N) 

1200  C  LOWER  AND  UPPER  BRANCH  ANALYSIS 
1210  10  IFiYP.GT . YLMAX(M) )  GO  TO  11 

1220  CALL  ROOTS! 1 f N> 41 »YP»A1UP>A2UP) 

1230  CALL  ROOTS ! 0 1 M » 41  *  YP » A1LP » A2LP ) 

1240  CALL  ROOTS! 1 *H»  41 » YH*A1UN*A2UM) 

1250  CALL  ROOTS! 0»M* 41 »YM»A1LH»A2LM) 

1260  ALB!M)*!AlUP4AlLP+AlUH+AlLH)/4. 

1270  BLB!M)=( AlLP-AlUM)/2. 

1280  ELB!M)-100.*DLB(M)/ALB!H) 

1290  AUB!M)=! A2LP+A2UP4A2LH+A2UM)/4  * 

1300  DUB ! N )  =  <  A2UH-A2LP ) / 2 . 

1310  EUB!M>=100.*DUB!M)/AUB!M> 

1320  11  CONTINUE 

1330  C  PRINT  RESULTS  AND  CONTINUE  BATA  PROCESSING  LOOP 
1340  WRITEI6.406) 

1350  WRITE! 6»  407 )  TR >L »C» D* Y » ( ALB! M) » AUB! M) » AOB(H) tM=l *2) 

1360  WRITE!6>407)  DTR»DL *DC»DD»DY* !DLB! H) >DUB!M ) » D0B!H) >M*1 ,2 ) 

1370  WRITE16.407)  ET » ELiEC>ED*EY» (ELB(M) >EUB(M) »EQB(H) >M=1 »2) 

1380  GO  TO  6 

1290  C 
HOC  C 

1410  100  FORMAT (5E10.0»9X» A 1 ) 

1420  101  F0RMAT18E10.0) 

1430  200  FORMAT ! 1H1 */4X»#A32C0DE  RESULTS*  /> 

1440  201  F0RMAT!8X>*NLIMITS*>UX»*N1=*>1PE10.3»3X>*N2=*>E10.3  ) 

1450  202  FORMAT !8X»*KLIMITS*r 11X**K1S*»1PE10.3>3X»*K2S**E10.3  /> 

1460  203  FORMAT! 8X»*WAVELENGTH*SfF8.4r*  MICRON*) 

1470  204  FORMAT (8X«*DISTRIBUTIQNs*« 11X»A1  ///) 

1480  300  FORMAT (4X»*Q/A32  BOUNDS  !UNIT=1/MU>*  /) 

1490  301  FORMAT (5X»*I*>4X»*A32!MU)*»7X»*X32*»7X»*AVE+SIG*>5X»*AVE-SIG*»7X. 

1500  X  *MAX*»9X»*M1N*  /) 

1510  302  FORMAT (I6*1P6E12.3) 

1520  400  FORMAT ! 1H1 >/3X»*A32  AND  TOTAL  ERROR  RESULTS*  /> 

1530  401  FORMAT ! 6X» *FIRST  ROW*VALUE  LB=LOWER  BRANCH*) 

1540  402  F0RMAT!6X»*SEC0ND  ROW=RMS  ERROR  UB*UPPER  BRANCH*) 

1550  403  FORMAT! 6X>*THIRD  ROW*PERCENT  ERROR  OB»UNOETERMINED  BRANCH*) 

1560  404  FORMAT! /6X»*ZER0*N0T  APPLICABLE*  /> 

1570  405  F0RMAT(5X**TRANS*»5X»*L!CM)*»4X»*C!G/CM3)*>2X»*D!G/CM3)*ilXr 

1580  X  *Q/A!MU-1)*»2X»*A32!MU)  RMS  ERROR  B0UNDS*»3X»*A32!MU)  FULL  ERROR 
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1590  XBOUNDS*  /56X»«tB*»7X»*UB*f7X>*0B*f7X»*lB*>7X»*UBt»?Xf.«0B*  /> 
1600  406  FQRHAT</> 

1610  407  F0RMAT(2X*1P5E10.3»6E9.2) 

1620  END 


1630  SUBROUTINE  REABQ(DISTR) 

1640  C 

1650  C  READ  QFILE  DATA 

1660  C 

1670  REAL  NINBEX<9) >KINDEX<6> 

1680  DIMENSION  X32<41) 

1690  DIMENSION  BBAR<9»6»41> »0SI6(9»6»41 ) »0MIN(9»6»41 > >QMAX(9»6r 41 ) 

1700  COMMON/BDATA/X32.NN»NKtNINDEX*KINBEX»BBAR*QSIQ*QMINfQHAX 

1710  C. 

1720  IF(DISTR.EQ.IHU)  M*2 

1730  IF(DISTR.EQ.IHB)  M*3 

1740  READ<M« 100)  < X32 < I ) » 1*1 »41 > 

1750  READ(M» 101 >  NNiNK 

1760  DO  1  N- 1 *NN 

1770  DO  1  K=1 *  NK 

1780  READ(M.IOO)  NINDEX(N) fKINBEX(K) 

1790  READ(M.IOO)  <0BAR(N,K.I>.I=1»41> 

1800  READ (H* 100)  <QSIG(N»K»I) .1*1.41) 

1810  READ(HtlOO)  <QMIN<N»K» I ) » 1*1 > 41 ) 

1820  READ(HflOO)  (RMAX(N,K» I ) » 1  =  1 »41 ) 

1830  1  CONTINUE 

1840  RETURN 

1850  ('• 

1860  100  FORMAT(BEIO.O) 

1870  101  FORMAT (21 10) 

I860  END 


SUBROUTINE  OINT(NINDEX. KINDEX) 


1890 
1900  f, 

1910  C  INTERPOLATE  ON  OFILE  DATA  TO  GET  QDATA  AT  NINDEX  AND  KINDEX 

1920  C 
1930  C 

1940  DIMENSION  FBAR<2) »FSIB<2> iFMIN(2> »FNAX(2) 

1950  REAL  NINDEX > KINDEX »NFIL£< 9) »KFILE (6) 

1960  DIMENSION  6BAR<9»6»41)»6SIG<9»6»41>tGMIN(9.6»41)»GMAX(9.6*41) 

1970  COMMON/QDATA/XXXXX(41 ) »NN»NK»NFILE*KFILE  *GBARr6SIGt GHINfGHAX 

1980  COMMON/QFUNC/OBAR( 41 )?0SI6(41)»QMIN(41> >QMAX(41 > 

1990  C 

2000  C  LOCATE  FILE  ENTRY  POINTS 

2010  DO  1  N=1.NN 

2020  IF(NFILE(N),GT. NINDEX)  GO  TO  2 

2030  1  CONTINUE 

2040  N=NN 

2050  2  DO  3  K=1 »NK 

2060  IF<KFILE(K).6T. KINDEX)  GO  TO  4 

2070  3  CONTINUE 

2080  K=NK 

2090  C  SIZE  PARAMETER  LOOP 

2100  4  DO  6  I  =  1 f  4 1 

2110  C  INTERPOLATE  ON  N  'LINEAR-LINEAR) 

2120  DO  5  L=1 » 2 

2130  M=K-2+L 

2140  DEL*(NINDEX-NFILE(N-l))/(NFILE<N)-NFILE(N-t>) 

2150  FBAR<L)=GBAR(N-1iMfI)+(GBAR(N»M»I)-GBAR<N-1.M»I))*DEL 

2160  FSIG(L)=GSIG(N-1>M>I)4<GSI6(NfMfI)-GSIG(N-1tMfI>)«DEL 

2170  FMIN(L)=GMIN(N-1»M»I)+(GMIN(N.M.I)-GMIN(N-1»M.I))*DEL 

2180  5  FMAX(L) -GMAX(N-1 »M»I )+(GHAX(N»H»I >-GMAX(N-l »M» I ) )*TiEL 

2190  C  INTERPOLATE  ON  K  (LOG-LINEAR) 

2200  DEL=(KINDEX-KFILE(K-1 ) )/(KFILE(K)-KFILE(K-l ) ) 

2210  QBAR<  I  )*FBAR(  1  >*<FBAR<  2  )/FBAR<  1 )  XtDEL 

2220  QSIG(I)=FSIG(1)*(FSIG(2)/FSIG(1))**DEL 

2230  QMIN(I)=FMIN(1)*<FHIN(2)/FHIN(1) XtDEL 

2240  6  QMAX< I )=FMAX( 1 X(FMAX(2)/FMAX< 1 ) XtDEL 

2250  RETURN 

2260  END 
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SUBROUTINE  CURE IT ( IBOUNDiM»N>DELX* YMAX > 


100 
110  C 

120  C  QUADRATIC  CURVE  FIT  TO  Q/A32  VS  A32  RESULTS 

130  C 

HO  C0NM0N/B0UND/XXXXX<41)*BU<2f41)»BL(2»41) 

ISO  CQMM0N/C0£FF/ACU(2»40) »BCU<2>40 ) »CCU( 2»  40) « ACL (2* 40) iBCL( 2 » 40) ? 

160  X  CCL(2»40> 

170  C 

180  C  LOOP  TO  OBTAIN  QUADRATIC  COEFFICIENTS  FOR  ALL  A  INTERVALS  AND 
190  C  DETERMINE  ABSOLUTE  MAXIMUM 
200  YMAX=0. 

210  NN=N-2 

220  DO  5  1=1, NN 

230  IF < I BOUND. EQ 4 1 )  60  TO  1 

240  Y1=BL<M« I ) 

250  Y2*BL(M» 1  +  1 ) 

260  Y3=BL(M» 1+2) 

270  GO  TO  2 

280  1  Y1=BU(M» I ) 

290  Y2=BU(M» 1  +  1 ) 

300  Y3=BU(M» 1+2) 

310  2  A=AL0G<Y1*Y3/Y2*t2)/(2.*DELX«2> 

320  B=AL0G<Y2**4/(Y3*Y1M3)>/<2.*DELX> 

330  C=ALOG( Y1 ) 

340  IFIA.GE.O.)  GO  TO  3 

350  TEST=-B/(2.*A) 

360  IFITEST.LT.  0.)  60  TO  3 

370  IF(TEST.GT.DELX)  60  TO  3 

380  YE=EXP<C-B*B/<4.*A>) 

390  IFtYE.GT . YMAX )  YMAX=YE 

400  3  CONTINUE 

410  IF < IBOUND.EQ. 1 )  GO  TO  4 

420  ACL(M,I)=A 

430  BCL ( H» I )=B 

440  CCL<M» I )=C 

450  GO  TO  5 

460  4  ACU(M, [ )  -A 

470  BCU(M» I )”B 

480  CCU(M* I )=C 

490  5  CONTINUE 

500  IF < IBOUND.EQ. 1 >  60  TO  6 

510  ACL(M»N-1 '=ACL'MiNN> 

520  BCUN,N-1>=BCL'M,NN> 

530  CCL ( M  >  N- 1 ) =CCL  <  M » NN ) 

540  RETURN 

550  6  ACU(M.N-1)=ACU<M»NN) 

560  BCU(M.N-1)«BCU(M*NN) 

570  CCU(M»N-l)=CCU(MiNN) 

580  RETURN 

590  END 
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SUBROUTINE  PLOTBI YUMAX) 
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PLOT  UPPER  AND  LONER  0/A32  VS  A32  CURVES 

DIMENSION  AP(41)»YU1<41 ) t YU2( 41 )»YL1(41)»YL2(41) 

DIMENSION  YUHAX<2> 

COMMON/BOUND/ A32 ( 4 1 > » BU ( 2 1 4 1 ) » BL  <  2 1 4 1 ) 

SET  UP  PLOT  LIMITS 
NP0INTS=4 l 

AMIN=IFIX( AL0G10I A32< 1 )>>—!* 

AMAX=IFIX(AL0G10(A32<41 ) )  >41, 

Y1MAX=IFIX( AL0610I YUMAXI 1 ) ) )  +  l . 

YlMIN=YlMAX-3. 

Y2MAX=IFIX(AL0G10(YUMAX(2)>)+1. 

Y2MIN=Y2MAX-3. 

SET  UP  PLOT  VARIABLES 
DO  1  I  =  1 1 4 i 
AP(  I )*AL0610( A32( I ) ) 

YUl(I)*AL0610<BUa*I)) 

YL1 ( I >=AL0G10<BL< 1 »I ) ) 

YU2(I)=AL0G10(BU(2»I)i 
YL2(I)=AL0G10(BL'2»I> ) 

IF(YU1<I).LT.Y1MIN)  YU1<I>*Y1MIN 
IF(YLI'I).LT.YIMIN)  YL1<I>*Y1MIN 
IF<YU2<I).LT.Y2MIN)  YU2(I>*Y2MIN 
IF(YL2(I).LT,Y2M1N)  YL2(I)*Y2MIN 
1  CONTINUE 

PLOT  YUl(I)  AND  YL1  ( I )  VS.  APU)  I«l>2* . . . >NPOINTS  ON  LINEAR- 
LINER  PLOT  TO  GET  UPPER  AND  LONER  BOUND  CURVES  APPROPRIATE  TO  RMS 
ERROR  IN  SIZE  DISTRIBUTION 

PLOT  YU2(I)  AND  YL2 < I >  VS.  AP(I)  TO  BET  CORRESPONDING  RESULT  FOR 
FULL  RANGE  ERROR  IN  SIZE  DISTRIBUTION 

SUGGESTED  X-AXIS  TITLE  —  LOGIO  A32<HICRON> 

SUGGESTED  Y-AXIS  TITLE  — -  LOGIO  0/A32I 1/H1CR0N > 

SUGGESTED  TITLE  FOR  FIRST  PLOT  —  ERROR  BOUNDS  FOR  RMS  SIZE  ERROR 
SUGGESTED  TITLE  FOR  SECOND  PLOT  —  ERROR  BOUNDS  FOR  FULL-RANGE  SIZE 

ERROR 

DELETE  ABOVE  COMMENTS  AND  INSERT  YOUR  PLOTTING  ROUTINES  HERE 
RETURN 
END 


2770  SUBROUTINE  ROOTS! IB0UND,M»N,Y,A1 ,A2> 

2780  C 

2790  C  BETERHINE  LOWEST  AND  HIGHEST  INTERSECTIONS  OF  Y  WITH 

2800  C  BOUND  CURVE 

2810  C 

2820  DIHENSION  ARUO) 

2820  C0MM0N/B0UND/A32  1  41 ) *  XXXXX  1 1 64 ) 

2840  CQHMQN/COEFF /ACU ( 2 » 40 ) *  BCU<  2  *  40  >  »CCU<  2 » 40 ) ,  ACL ( 2»  40 ) » BCL ( 2»  40 ) 

2850  X  CCL(2r 40) 

2860  DATA  DELX/O . 230258510/ 

2870  C 

2880  C  SEARCH  FOR  ROOTS  IN  ALL  A32  INTERVALS 


2890 

NR=0 

2900 

NN=N-1 

2910 

DO  4  1=1, NN 

2920 

IF< IB0UND.E9. 1 )  GO  TO  1 

2930 

A=ACL( H» I > 

2940 

B=BCL1N» I > 

2950 

OCCL1M,  I > 

2960 

GO  TO  2 

2970 

1 

A=ACU(M , I ) 

2980 

B=BCU(W, I ) 

2990 

C=CCU(H, I ) 

3000 

2 

TEST=B*B-4.*A*<C-AL0G1Y)) 

3010 

IF1TEST.LT.0.)  GO  TO  4 

3020 

TEST1=-<B-SGRT<TEST))/<2.*A) 

3030 

TEST2=-1B+SGRT <  TEST) )/<2.*A) 

3040 

IFITESTl.LT. 0.)  GO  TO  3 

3050 

IF1TEST1 .GE.DELX)  GO  TO  3 

3060 

NR=NR+1 

3070 

AR<NR)=A32(I)*EXP1TEST1) 

3090 

3 

IF1TEST2.LT.0.)  GO  TO  4 

3090 

IF1TEST2. GE.DELX)  GO  TO  4 

3100 

NR=NR+1 

3110 

AR(NR)=A32(I'tEXP(TEST2) 

3120 

4 

CONTINUE 

2130  C  FIND  LOWEST  AND  HIGHEST  ROOTS 
3140  IF(NR.NE.O)  GO  TO  5 

3150  A1=0. 

3160  A2=0. 

3130  RETURN 

2180  5  IF1NR.NE.1;  GO  TO  6 

3190  A1=0 . 

3200  A2=AR< 1 ) 

3210  RETURN 

3220  6  A 1 = 1 .E99 

3230  A2=0. 

3240  DO  7  L=1 »NR 

3250  IFIARID.LT.AI)  A1=AR(L) 
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3240  IF(AR(L) «6T. A2>  A2*AR(L) 

3270  7  CONTINUE 

32HO  RETURN 

3290  END 
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END 
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